{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Regenerating Data from\n",
    "# [J.T. Gostick et al. / JPS 173 (2007) 277–290](http://www.sciencedirect.com/science/article/pii/S0378775307009056)\n",
    "\n",
    "## Getting Started\n",
    "\n",
    "In this tutorial, we will regenerate data from J.T. Gostick's 2007 paper [[1]](http://www.sciencedirect.com/science/article/pii/S0378775307009056). This will both show that OpenPNM can recreate results accurately, and will also show some more specific uses of OpenPNM. While this paper deals with both SGL and Toray GDLs, we will deal only with SGL.\n",
    "\n",
    "There will be a general layout to complete this simulation:\n",
    "\n",
    "1. Set up network\n",
    "2. Set up geometry and geometrical methods\n",
    "3. constrict throat's by a constriction factor\n",
    "4. Set up phases and methods\n",
    "5. Set up phase physics and methods\n",
    "6. Run invasion percolation\n",
    "7. Run Stokes and Fickian algorithms\n",
    "8. generate effective permeability and effective diffusivity values at different saturations\n",
    "9. plot generated data\n",
    "\n",
    "We first import the openpnm code and some other useful modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:14.735456Z",
     "iopub.status.busy": "2021-06-24T11:24:14.733733Z",
     "iopub.status.idle": "2021-06-24T11:24:15.501340Z",
     "shell.execute_reply": "2021-06-24T11:24:15.499819Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import openpnm as op\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "import matplotlib.pyplot as plt\n",
    "import openpnm.models as mods\n",
    "%matplotlib inline\n",
    "np.random.seed(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up Network and Geometry\n",
    "\n",
    "To begin our simulation, we must first generate our SGL network and geometry.  This includes:\n",
    "\n",
    "1. creating a cubic network object and an SGL10 geometry object\n",
    "2. sending our geometry object our internal pores\n",
    "3. calculating values for throat and pore properties for both internal and boundary pores\n",
    "4. accounting for pores and throats that are too big (making maximum pore size the lattice parameter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:15.562656Z",
     "iopub.status.busy": "2021-06-24T11:24:15.561138Z",
     "iopub.status.idle": "2021-06-24T11:24:15.832415Z",
     "shell.execute_reply": "2021-06-24T11:24:15.833577Z"
    }
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "module 'openpnm.geometry' has no attribute 'Boundary'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_722488/1704912735.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     47\u001b[0m \u001b[0mPs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msgl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpores\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'*boundary'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     48\u001b[0m \u001b[0mTs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msgl\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfind_neighbor_throats\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpores\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mPs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'or'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 49\u001b[0;31m \u001b[0mboun\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mop\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgeometry\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBoundary\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnetwork\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msgl\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpores\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mPs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mthroats\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mTs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'boun'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m: module 'openpnm.geometry' has no attribute 'Boundary'"
     ]
    }
   ],
   "source": [
    "Lc = 40.5e-6\n",
    "# 1. Set up network\n",
    "sgl = op.network.Cubic(shape=[26, 26, 10], spacing=Lc, name='SGL10BA')\n",
    "sgl.add_boundary_pores()\n",
    "proj = sgl.project\n",
    "wrk = op.Workspace()\n",
    "wrk.settings['loglevel'] = 50\n",
    "# 2. Set up geometries\n",
    "Ps = sgl.pores('*boundary', mode='not')\n",
    "Ts = sgl.find_neighbor_throats(pores=Ps, mode='xnor', flatten=True)\n",
    "geo = op.geometry.GenericGeometry(network=sgl,pores=Ps,throats=Ts,name='geo')\n",
    "geo.add_model(propname='pore.seed',\n",
    "              model=mods.misc.random,\n",
    "              element='pore',\n",
    "              num_range=[0, 0.8834],\n",
    "              seed=None)\n",
    "geo.add_model(propname='throat.seed',\n",
    "              model=mods.misc.from_neighbor_pores,\n",
    "              prop='pore.seed',\n",
    "              mode='min')\n",
    "geo.add_model(propname='pore.diameter',\n",
    "              model=mods.geometry.pore_size.weibull,\n",
    "              shape=3.07,\n",
    "              loc=19.97e-6,\n",
    "              scale=1.6e-5)\n",
    "geo.add_model(propname='throat.diameter',\n",
    "              model=mods.geometry.throat_size.weibull,\n",
    "              shape=3.07,\n",
    "              loc=19.97e-6,\n",
    "              scale=1.6e-5)\n",
    "geo.add_model(propname='pore.area',\n",
    "              model=mods.geometry.pore_cross_sectional_area.sphere)\n",
    "geo.add_model(propname='pore.volume',\n",
    "              model=mods.geometry.pore_volume.sphere)\n",
    "geo.add_model(propname='throat.length',\n",
    "              model=mods.geometry.throat_length.ctc)\n",
    "geo.add_model(propname='throat.volume',\n",
    "              model=mods.geometry.throat_volume.cylinder)\n",
    "geo.add_model(propname='throat.area',\n",
    "              model=mods.geometry.throat_cross_sectional_area.cylinder)\n",
    "geo.add_model(propname='throat.surface_area',\n",
    "              model=mods.geometry.throat_surface_area.cylinder)\n",
    "geo.add_model(propname='throat.endpoints',\n",
    "              model=mods.geometry.throat_endpoints.spherical_pores)\n",
    "geo.add_model(propname='throat.conduit_lengths',\n",
    "              model=mods.geometry.throat_length.conduit_lengths)\n",
    "Ps = sgl.pores('*boundary')\n",
    "Ts = sgl.find_neighbor_throats(pores=Ps, mode='or')\n",
    "boun = op.geometry.Boundary(network=sgl, pores=Ps, throats=Ts, name='boun')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we move on to setting up our fluid and physics objects, we must constrict throats in the z and y direction by a factor (Gostick et al included this tightening of throats in only these two directions to create realistic anisotropy in the model).  For his SGL simulation, Gostick uses a constriction factor of .95.  Finally, because we have changed values for pore and throat diameters (first by accounting for pores and throats that are too big, and the finally constricting throats in the y and z directions), we must recalculate all pore and throat values relying on these diameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:15.841566Z",
     "iopub.status.busy": "2021-06-24T11:24:15.840142Z",
     "iopub.status.idle": "2021-06-24T11:24:16.091770Z",
     "shell.execute_reply": "2021-06-24T11:24:16.089244Z"
    }
   },
   "outputs": [],
   "source": [
    "throats = geo.throats()\n",
    "connected_pores = sgl.find_connected_pores(throats)\n",
    "x1 = [sgl['pore.coords'][pair[0]][0] for pair in connected_pores]\n",
    "x2 = [sgl['pore.coords'][pair[1]][0] for pair in connected_pores]\n",
    "same_x = [x - y == 0 for x, y in zip(x1,x2)]\n",
    "factor = [s*.95 + (not s)*1 for s in same_x]\n",
    "throat_diameters = sgl['throat.diameter'][throats]*factor\n",
    "geo['throat.diameter'] = throat_diameters\n",
    "geo.regenerate_models(exclude=['throat.diameter'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenPNM makes it very easy to visualize the network we have generated through the \"Visualization\" methods.  We can create vtk files to be viewed using ParaView (downloadable at http://www.paraview.org/download/ ). If we visualize our pore network model it would appear like this (the pores have been visualized using boxes- darker boxes are larger. Because the network is so big, visualization of the throats has been left out for clarity):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:16.099778Z",
     "iopub.status.busy": "2021-06-24T11:24:16.098279Z",
     "iopub.status.idle": "2021-06-24T11:24:16.775007Z",
     "shell.execute_reply": "2021-06-24T11:24:16.773637Z"
    }
   },
   "outputs": [],
   "source": [
    "import openpnm.io.VTK as iovtk\n",
    "iovtk.export_data(network=sgl, filename='network_SGL')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example is seen here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"http://i.imgur.com/fPZ8lZK.png\" style=\"width: 60%\" align=\"left\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the Phases and Physics\n",
    "Now we are ready to set up our phases (water and air) and the physics corresponding to each of these phases. openpnm has built in air and water phases, so we can use those. However, Gostick specifies using a water pore contact angle of 100, so we will reset this value after regenerating our fluids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:16.782581Z",
     "iopub.status.busy": "2021-06-24T11:24:16.781186Z",
     "iopub.status.idle": "2021-06-24T11:24:16.806010Z",
     "shell.execute_reply": "2021-06-24T11:24:16.807224Z"
    }
   },
   "outputs": [],
   "source": [
    "air = op.phases.Air(network = sgl, name = 'air')\n",
    "water = op.phases.Water(network = sgl, name = 'water')\n",
    "# Reset pore contact angle\n",
    "water['pore.contact_angle'] = 100.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to establish physical properties for our fluid objects. To do this, we will: 1) create physics objects associated with our fluids (by using StandardPhyics we don't have to add methods for calculating each property because they are already included) 2) use our regenerate_physics() method to calculate these properties. One physics object is required for each combination of phase and geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:16.815886Z",
     "iopub.status.busy": "2021-06-24T11:24:16.814389Z",
     "iopub.status.idle": "2021-06-24T11:24:17.199821Z",
     "shell.execute_reply": "2021-06-24T11:24:17.198565Z"
    }
   },
   "outputs": [],
   "source": [
    "phys_water = op.physics.Standard(network=sgl, phase=water, geometry=geo)\n",
    "phys_air = op.physics.Standard(network=sgl, phase=air, geometry=geo)\n",
    "phys_water_b = op.physics.Standard(network=sgl, phase=water, geometry=boun)\n",
    "phys_air_b = op.physics.Standard(network=sgl, phase=air, geometry=boun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running Ordinary Percolation, Fickian Diffusion, and Stokes Flow\n",
    "Gostick uses ordinary percolation to spread water through his GDL before calculating relative permeability and relative diffusivity.  This way, a graph showing the relationship between saturation and relative permeability and between saturation and relative diffusivity can be created.\n",
    "\n",
    "To run our ordinary percolation, we will:\n",
    "\n",
    "1. pick inlet pores\n",
    "2. create an Ordinary Percolation algorithm object\n",
    "3. setup our algorithm object\n",
    "4. run our algorithm object\n",
    "5. call results() so that occupancy of pores and throats for each fluid will can be set and multiphysics updated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:17.207818Z",
     "iopub.status.busy": "2021-06-24T11:24:17.206407Z",
     "iopub.status.idle": "2021-06-24T11:24:17.387500Z",
     "shell.execute_reply": "2021-06-24T11:24:17.386167Z"
    }
   },
   "outputs": [],
   "source": [
    "inlets = sgl.pores('bottom_boundary')\n",
    "used_inlets = [inlets[x] for x in range(0, len(inlets), 2)]\n",
    "\n",
    "OP_1 = op.algorithms.OrdinaryPercolation(project=proj)\n",
    "OP_1.set_inlets(pores=used_inlets)\n",
    "OP_1.settings.update({'phase': water.name,\n",
    "                      'pore_volume': 'pore.volume', \n",
    "                      'throat_volume': 'throat.volume'})\n",
    "OP_1.run(points=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This algorithm performed a start to finish simulation, which fully flooded the network. The 'results()' command can be used to update the phase occupancy values throughout the network. To save some computation, we will filter the invasion points so that relative transport properties can be calculated approximately every 5% increment in saturation. The OrdinaryPercolation object has a method to return the intrusion data as a named tuple of Capillary Pressure (Pcap) and Saturation of the non-wetting phase (Snwp)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:17.396424Z",
     "iopub.status.busy": "2021-06-24T11:24:17.394526Z",
     "iopub.status.idle": "2021-06-24T11:24:17.409158Z",
     "shell.execute_reply": "2021-06-24T11:24:17.410233Z"
    }
   },
   "outputs": [],
   "source": [
    "data = OP_1.get_intrusion_data()\n",
    "# Filter for evenly spaced sat inc. first and last\n",
    "filter_pc = [data.Pcap[0]]\n",
    "sat = [data.Snwp[0]]\n",
    "for i, pc in enumerate(data.Pcap):\n",
    "    if  data.Snwp[i] - sat[-1] > 0.05:\n",
    "        filter_pc.append(pc)\n",
    "        sat.append(data.Snwp[i])\n",
    "filter_pc.append(data.Pcap[-1])\n",
    "sat.append(data.Snwp[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define a helper function to update the phases and properties with the results of the OP algorithm. The multiphase conduit conductance model looks at the phase occupancy in the conduits made by the 1/2 pore - throat - 1/2 pore neighbor elements. When the mode is 'strict' the phase must occupy all three elements for the conduit to be considered open to flow for that phase. If the phase is not present in at least one of the elements in the conduit then the throat conductance is divided by 6 orders of magnitude. In this way the conductivity is severely reduced by the presence of the other phase and flow must go around, thus decreasing the permeability/diffusivity of the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:17.420382Z",
     "iopub.status.busy": "2021-06-24T11:24:17.418814Z",
     "iopub.status.idle": "2021-06-24T11:24:17.421862Z",
     "shell.execute_reply": "2021-06-24T11:24:17.422901Z"
    }
   },
   "outputs": [],
   "source": [
    "def update_phase_and_phys(results):\n",
    "    water['pore.occupancy'] = results['pore.occupancy']\n",
    "    air['pore.occupancy'] = 1-results['pore.occupancy']\n",
    "    water['throat.occupancy'] = results['throat.occupancy']\n",
    "    air['throat.occupancy'] = 1-results['throat.occupancy']\n",
    "    # Add multiphase conductances\n",
    "    mode='strict'\n",
    "    phys_air.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                       propname='throat.conduit_diffusive_conductance',\n",
    "                       throat_conductance='throat.diffusive_conductance',\n",
    "                       mode=mode)\n",
    "    phys_water.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                         propname='throat.conduit_diffusive_conductance',\n",
    "                         throat_conductance='throat.diffusive_conductance',\n",
    "                         mode=mode)\n",
    "    phys_air.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                       propname='throat.conduit_hydraulic_conductance',\n",
    "                       throat_conductance='throat.hydraulic_conductance',\n",
    "                       mode=mode)\n",
    "    phys_water.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                         propname='throat.conduit_hydraulic_conductance',\n",
    "                         throat_conductance='throat.hydraulic_conductance',\n",
    "                         mode=mode)\n",
    "    phys_air_b.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                         propname='throat.conduit_diffusive_conductance',\n",
    "                         throat_conductance='throat.diffusive_conductance',\n",
    "                         mode=mode)\n",
    "    phys_water_b.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                           propname='throat.conduit_diffusive_conductance',\n",
    "                           throat_conductance='throat.diffusive_conductance',\n",
    "                           mode=mode)\n",
    "    phys_air_b.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                         propname='throat.conduit_hydraulic_conductance',\n",
    "                         throat_conductance='throat.hydraulic_conductance',\n",
    "                         mode=mode)\n",
    "    phys_water_b.add_model(model=mods.physics.multiphase.conduit_conductance,\n",
    "                           propname='throat.conduit_hydraulic_conductance',\n",
    "                           throat_conductance='throat.hydraulic_conductance',\n",
    "                           mode=mode)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following call will get the pore and throat phase occupancy which is an array of 1s and 0s representing that the phase occupies a particular pore or throat, update the phase objects and and multiphase conductanct models to the physics objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:17.429440Z",
     "iopub.status.busy": "2021-06-24T11:24:17.428205Z",
     "iopub.status.idle": "2021-06-24T11:24:17.461292Z",
     "shell.execute_reply": "2021-06-24T11:24:17.459966Z"
    }
   },
   "outputs": [],
   "source": [
    "update_phase_and_phys(OP_1.results(Pc=1e3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step will be to calculate effective diffusivity and permeability at different saturations.  Note that we want to run Fickian diffusion and Stokes flow algorithms at different points within our ordinary percolation process.\n",
    "\n",
    "The rest of our code will exist within a loop updating our network to different stages of percolation, so that we may view our relative diffusivity and permeability at different points of saturation.\n",
    "\n",
    "Before we add in the loop aspect, we will walk through the code that will be inside the loop.\n",
    "\n",
    "Note that we want the algorithms that are single phase (where only the specified fluid exists in the network) to help us make our permeability and diffusivity values relative.  Any algorithm that is single phase will use the hydraulic or diffusive conductances before we recalculated based on occupancy.  This calls for our conductance parameter to be 'hydraulic_conductance' or 'diffusive_conductance' instead of 'conduit_hydraulic_conductance' or 'conduit_diffusive_conductance'.\n",
    "\n",
    "The need for all these different algorithms can be made clearer by the equation relating effective permeability to the absolute permeability and relative permeability:\n",
    "\n",
    "$$K_{eff, p}(s_p) = K*K_{r, p}(s_p)$$\n",
    "\n",
    "|Symbol|Description|\n",
    "|-------------------|-------------------------------------------------------------|\n",
    "|$$K_{eff, p}(s_p)$$|effective permeability of phase p as a function of saturation|\n",
    "|$$K$$              |absoulte permeability (or single phase permeability)         |\n",
    "|$$K_{r, p}(s_p)$$  |relative permeability of phase p as a function of saturation |\n",
    "\n",
    "Therefore, relative permeability can be found by dividing the effective permeability by the absolute permeability.  Thus the need for a single phase algorithm (absolute permeability) for every multi phase algorithm (effective permeability).\n",
    "\n",
    "The same goes for relative diffusivity, which has an very similar equation that looks like this:\n",
    "$$D_{eff, p}(s_p) = D*D_{r, p}(s_p)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:17.479626Z",
     "iopub.status.busy": "2021-06-24T11:24:17.477977Z",
     "iopub.status.idle": "2021-06-24T11:24:19.445664Z",
     "shell.execute_reply": "2021-06-24T11:24:19.446188Z"
    }
   },
   "outputs": [],
   "source": [
    "perm_air = {'0': [], '1': [], '2': []}\n",
    "diff_air = {'0': [], '1': [], '2': []}\n",
    "perm_water = {'0': [], '1': [], '2': []}\n",
    "diff_water = {'0': [], '1': [], '2': []}\n",
    "\n",
    "max_Pc = max(OP_1['throat.invasion_pressure'])\n",
    "\n",
    "num_seq = 20\n",
    "pore_volumes = sgl['pore.volume']\n",
    "throat_volumes = sgl['throat.volume']\n",
    "totV = np.sum(pore_volumes) + np.sum(throat_volumes)\n",
    "\n",
    "K_air_single_phase = [None, None, None]\n",
    "D_air_single_phase = [None, None, None]\n",
    "K_water_single_phase = [None, None, None]\n",
    "D_water_single_phase = [None, None, None]\n",
    "bounds = [['left', 'right'], ['front', 'back'], ['top', 'bottom']]\n",
    "\n",
    "for bound_increment in range(len(bounds)):\n",
    "    # Run Single phase algs effective properties\n",
    "    BC1_pores = sgl.pores(labels=bounds[bound_increment][0]+'_boundary')\n",
    "    BC2_pores = sgl.pores(labels=bounds[bound_increment][1]+'_boundary')\n",
    "    \n",
    "    # Effective permeability : air\n",
    "    sf_air = op.algorithms.StokesFlow(network=sgl, phase=air)\n",
    "    sf_air.settings['conductance'] = 'throat.hydraulic_conductance'\n",
    "    sf_air.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "    sf_air.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "    sf_air.run()\n",
    "    K_air_single_phase[bound_increment] = sf_air.calc_effective_permeability()\n",
    "    proj.purge_object(obj=sf_air)\n",
    "    \n",
    "    # Effective diffusivity : air\n",
    "    fd_air = op.algorithms.FickianDiffusion(network=sgl,phase=air)\n",
    "    fd_air.settings['conductance'] = 'throat.diffusive_conductance'\n",
    "    fd_air.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "    fd_air.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "    fd_air.run()\n",
    "    D_air_single_phase[bound_increment] = fd_air.calc_effective_diffusivity()\n",
    "    proj.purge_object(obj=fd_air)\n",
    "    \n",
    "    # Effective permeability : water\n",
    "    sf_water = op.algorithms.StokesFlow(network=sgl, phase=water)\n",
    "    sf_water.settings['conductance'] = 'throat.hydraulic_conductance'\n",
    "    sf_water.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "    sf_water.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "    sf_water.run()\n",
    "    K_water_single_phase[bound_increment] = sf_water.calc_effective_permeability()\n",
    "    proj.purge_object(obj=sf_water)\n",
    "    \n",
    "    # Effective diffusivity : water\n",
    "    fd_water = op.algorithms.FickianDiffusion(network=sgl,phase=water)\n",
    "    fd_water.settings['conductance'] = 'throat.diffusive_conductance'\n",
    "    fd_water.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "    fd_water.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "    fd_water.run()\n",
    "    D_water_single_phase[bound_increment] = fd_water.calc_effective_diffusivity()\n",
    "    proj.purge_object(obj=fd_water)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can repeat the algorithms at each filtered pressure. This process takes about 1 minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:19.456292Z",
     "iopub.status.busy": "2021-06-24T11:24:19.455691Z",
     "iopub.status.idle": "2021-06-24T11:24:47.395709Z",
     "shell.execute_reply": "2021-06-24T11:24:47.395129Z"
    }
   },
   "outputs": [],
   "source": [
    "for Pc in filter_pc:\n",
    "    update_phase_and_phys(OP_1.results(Pc=Pc))\n",
    "    print('-' * 80)\n",
    "    print('Pc', Pc)\n",
    "    for bound_increment in range(len(bounds)):\n",
    "        BC1_pores = sgl.pores(labels=bounds[bound_increment][0]+'_boundary')\n",
    "        BC2_pores = sgl.pores(labels=bounds[bound_increment][1]+'_boundary')\n",
    "\n",
    "        # Multiphase\n",
    "        sf_air = op.algorithms.StokesFlow(network=sgl,phase=air)\n",
    "        sf_air.settings['conductance'] ='throat.conduit_hydraulic_conductance'\n",
    "        sf_water = op.algorithms.StokesFlow(network=sgl,phase=water)\n",
    "        sf_water.settings['conductance'] = 'throat.conduit_hydraulic_conductance'\n",
    "\n",
    "        fd_air = op.algorithms.FickianDiffusion(network=sgl,phase=air)\n",
    "        fd_air.settings['conductance'] = 'throat.conduit_diffusive_conductance'\n",
    "        fd_water = op.algorithms.FickianDiffusion(network=sgl,phase=water)\n",
    "        fd_water.settings['conductance'] = 'throat.conduit_diffusive_conductance'\n",
    "\n",
    "        #BC1\n",
    "        sf_air.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "        sf_water.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "        fd_air.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "        fd_water.set_value_BC(values=0.6, pores=BC1_pores)\n",
    "\n",
    "        #BC2\n",
    "        sf_air.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "        sf_water.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "        fd_air.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "        fd_water.set_value_BC(values=0.2, pores=BC2_pores)\n",
    "\n",
    "        # Run Multiphase algs\n",
    "        sf_air.run()\n",
    "        sf_water.run()\n",
    "        fd_air.run()\n",
    "        fd_water.run()\n",
    "\n",
    "        Keff_air_mphase = sf_air.calc_effective_permeability()\n",
    "        Deff_air_mphase = fd_air.calc_effective_diffusivity()\n",
    "        Keff_water_mphase = sf_air.calc_effective_permeability()\n",
    "        Deff_water_mphase = fd_water.calc_effective_diffusivity()\n",
    "\n",
    "        Kr_eff_air = Keff_air_mphase / K_air_single_phase[bound_increment]\n",
    "        Kr_eff_water = Keff_water_mphase / K_water_single_phase[bound_increment]\n",
    "        Dr_eff_air = Deff_air_mphase / D_air_single_phase[bound_increment]\n",
    "        Dr_eff_water = Deff_water_mphase / D_water_single_phase[bound_increment]\n",
    "\n",
    "        perm_air[str(bound_increment)].append(Kr_eff_air)\n",
    "        diff_air[str(bound_increment)].append(Dr_eff_air)\n",
    "        perm_water[str(bound_increment)].append(Kr_eff_water)\n",
    "        diff_water[str(bound_increment)].append(Dr_eff_water)\n",
    "        \n",
    "        proj.purge_object(obj=sf_air)\n",
    "        proj.purge_object(obj=sf_water)\n",
    "        proj.purge_object(obj=fd_air)\n",
    "        proj.purge_object(obj=fd_water)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot the results including those from the paper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:24:47.439109Z",
     "iopub.status.busy": "2021-06-24T11:24:47.423805Z",
     "iopub.status.idle": "2021-06-24T11:24:47.776034Z",
     "shell.execute_reply": "2021-06-24T11:24:47.777163Z"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib.font_manager import FontProperties\n",
    "%matplotlib inline\n",
    "\n",
    "# Data points taken directly from Gostick's graphs using GraphClick\n",
    "gostick_saturation_1 = [0.008, 0.04, 0.093, 0.14, 0.193, 0.246, 0.293, 0.337, 0.395, 0.442, 0.496,\n",
    "                        0.542, 0.59, 0.641, 0.687, 0.748, 0.793, 0.838, 0.894, 0.945, 0.986]\n",
    "gostick_perm_air_case1 = [0.917, 0.821, 0.68, 0.568, 0.466, 0.366, 0.286, 0.204, 0.144, 0.096, 0.051, 0.024,\n",
    "                          0.003, -1.08E-04, -1.96E-04, -3.12E-04, -3.97E-04, -4.84E-04, -5.90E-04, 0.002, 0.002]\n",
    "gostick_saturation_2 = [0.99, 0.899, 0.847, 0.802, 0.75, 0.701, 0.645, 0.594, 0.546, 0.497, 0.449,\n",
    "                        0.398, 0.348, 0.298, 0.245, 0.196, 0.147, 0.094, 0.044, 0.003]\n",
    "gostick_perm_water = [0.935, 0.774, 0.709, 0.664, 0.618, 0.572, 0.514, 0.461, 0.401, 0.347,\n",
    "                        0.284, 0.211, 0.145, 0.084, 0.044, 0.024, 0.012, 0.001, 0.001, 0.001]\n",
    "\n",
    "gostick_saturation_3 =[0.006, 0.05, 0.102, 0.151, 0.199, 0.247, 0.297, 0.348, 0.399, 0.447, 0.496,\n",
    "                    0.546, 0.597, 0.645, 0.699, 0.75, 0.798, 0.846, 0.899, 0.949, 0.983]\n",
    "gostick_diff_air_case1 = [0.939, 0.836, 0.725, 0.626, 0.531, 0.442, 0.353, 0.27, 0.203, 0.14, 0.085, 0.048,\n",
    "                          0.008, 5.49E-04, 4.48E-04, 3.50E-04, 2.59E-04, 1.67E-04, 0.003, 0.003, 0.003]\n",
    "gostick_saturation_4 = [0.985, 0.946, 0.898, 0.846, 0.795, 0.749, 0.695, 0.643, 0.596, 0.545, 0.496, 0.448,\n",
    "                        0.396, 0.346, 0.298, 0.251, 0.196, 0.146, 0.094]\n",
    "gostick_diff_water = [0.941, 0.901, 0.853, 0.809, 0.756, 0.7, 0.638, 0.569, 0.503, 0.428, 0.36, 0.291, 0.214, 1.48E-01,\n",
    "                      8.00E-02, 4.50E-02, 2.30E-02, 1.60E-02, 0.005]\n",
    "\n",
    "fontP = FontProperties()\n",
    "fontP.set_size('small')\n",
    "# Setting up subplots\n",
    "fig = plt.figure(figsize=(6, 10), dpi=80, facecolor='w', edgecolor='k')\n",
    "ax1 = fig.add_subplot(211)   #top\n",
    "ax2 = fig.add_subplot(212)   #bottom\n",
    "\n",
    "x_values1 = [x/20 for x in range(21)]\n",
    "z = '.75'\n",
    "\n",
    "\n",
    "# Plots for subplot1 - strict permeability\n",
    "p1, = ax1.plot(sat, perm_water['0'], color = 'k', linestyle = '-', marker = 'o')\n",
    "p2, = ax1.plot(sat, perm_water['1'], color = z, linestyle = '-', marker = 'o')\n",
    "p3, = ax1.plot(sat, perm_water['2'], color = 'b', linestyle = '-', marker = 'o')\n",
    "p4, = ax1.plot(sat, perm_air['0'], color = 'k', linestyle = '-', marker = '^')\n",
    "p5, = ax1.plot(sat, perm_air['1'], color = z, linestyle = '-', marker = '^')\n",
    "p6, = ax1.plot(sat, perm_air['2'], color = 'b', linestyle = '-', marker = '^')\n",
    "p10, = ax1.plot(x_values1, [x**(3) for x in x_values1], 'k--')\n",
    "ax1.plot(x_values1, [(1-x)**(3) for x in x_values1], 'k--')\n",
    "gs1, = ax1.plot(gostick_saturation_1, gostick_perm_air_case1, color = 'r', linestyle = '-', marker = 'D')\n",
    "gs2, = ax1.plot(gostick_saturation_2, gostick_perm_water, color = 'r', linestyle = '-', marker = 'o')\n",
    "ax1.set_ylabel('permeability')\n",
    "ax1.set_xlabel(\"saturation\")\n",
    "ax1.set_ylim([0,1])\n",
    "ax1.set_xlim([0,1])\n",
    "\n",
    "# Need to work on legend to match up with the right things\n",
    "lgd1 = ax1.legend([p1, p2, p3, p4, p5, p6, p10, gs1, gs2],\n",
    "                  [\"KrWater,x\", \"KrWater,y\", \"KrWater,z\",\n",
    "                   \"KrAir,x\",\"KrAir,y\",\"KrAir,z\", \"a = 3\", \n",
    "                   \"Gostick et al \\n KrAir,x (case 1)\", \n",
    "                   \"Gostick et al \\n KrWater,x\"], \n",
    "                  loc='center left', bbox_to_anchor=(1, 0.5), prop=fontP)\n",
    "\n",
    "# Plots for subplot4 - diffusivity\n",
    "p11, = ax2.plot(sat, diff_water['0'], color = 'k', linestyle = '-', marker = 'o')\n",
    "p12, = ax2.plot(sat, diff_water['1'], color = z, linestyle = '-', marker = 'o')\n",
    "p13, = ax2.plot(sat, diff_water['2'], color = 'b', linestyle = '-', marker = 'o')\n",
    "p14, = ax2.plot(sat, diff_air['0'], color = 'k', linestyle = '-', marker = '^')\n",
    "p15, = ax2.plot(sat, diff_air['1'], color = z, linestyle = '-', marker = '^')\n",
    "p16, = ax2.plot(sat, diff_air['2'], color = 'b', linestyle = '-', marker = '^')\n",
    "p20, = ax2.plot(x_values1, [x**(2) for x in x_values1], 'k--')\n",
    "ax2.plot(x_values1, [(1-x)**(2) for x in x_values1], 'k--')\n",
    "gs3, = ax2.plot(gostick_saturation_3, gostick_diff_air_case1, color = 'r', linestyle = '-', marker = 'D')\n",
    "gs4, = ax2.plot(gostick_saturation_4, gostick_diff_water, color = 'r', linestyle = '-', marker = 'o')\n",
    "ax2.set_ylabel('diffusivity')\n",
    "ax2.set_xlabel(\"saturation\")\n",
    "ax2.set_ylim([0,1])\n",
    "ax2.set_xlim([0,1])\n",
    "\n",
    "lgd2 = ax2.legend([p11, p12, p13, p14, p15, p16, p20, gs3, gs4],\n",
    "                  [\"DrWater,x\", \"DrWater,y\", \"DrWater,z\",\n",
    "                   \"DrAir,x\",\"DrAir,y\",\"DrAir,z\", \"a = 2\", \n",
    "                   \"Gostick et al \\n DrAir,x (case 1)\", \n",
    "                   \"Gostick et al \\n DrWater,x\"], \n",
    "                  loc='center left', bbox_to_anchor=(1, 0.5), prop=fontP)\n",
    "\n",
    "fig.subplots_adjust(left=0.13, right=.7, top=0.95, bottom=0.05)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrepancies with Gostick's simulation\n",
    "Several things contribute to slight differences between this simulation and that produced by Gostick et al in their 2007 paper.  These include:\n",
    "\n",
    "1. Lack of pore size correlation\n",
    "2. Lack of late pore filling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Acknowledgements\n",
    "The OpenPNM team would like to thank Jackie Lunger (Materials Science and Engineering, University of Toronto, 1T7) for her excellent work in developing this example."
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
